{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tue Jul 23 11:31:36 PDT 2019\r\n"
     ]
    }
   ],
   "source": [
    "import numpy as np, matplotlib.pyplot as plt, pandas as pd\n",
    "pd.set_option('display.max_rows', 10)\n",
    "!date\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Consistent models in DisMod-MR from Vivarium artifact draw\n",
    "\n",
    "Take i, r, f, p from a Vivarium artifact, and make a consistent version of them.  See how it compares to the original."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(123456)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# if dismod_mr is not installed, it should possible to use\n",
    "# !conda install --yes pymc\n",
    "# !pip install dismod_mr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import dismod_mr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# you also need one more pip installable package\n",
    "# !pip install vivarium_public_health"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import vivarium_public_health"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Consistent fit with all data\n",
    "\n",
    "Let's start with a consistent fit of the simulated PD data.  This includes data on prevalence, incidence, and SMR, and the assumption that remission rate is zero.  All together this counts as four different data types in the DisMod-II accounting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from vivarium_public_health.dataset_manager import Artifact\n",
    "art = Artifact('/share/costeffectiveness/artifacts/obesity/obesity.hdf')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[EntityKey(metadata.keyspace),\n",
       " EntityKey(metadata.versions),\n",
       " EntityKey(metadata.locations),\n",
       " EntityKey(population.demographic_dimensions),\n",
       " EntityKey(population.structure),\n",
       " EntityKey(cause.all_causes.cause_specific_mortality),\n",
       " EntityKey(population.theoretical_minimum_risk_life_expectancy),\n",
       " EntityKey(cause.ischemic_heart_disease.restrictions),\n",
       " EntityKey(cause.ischemic_stroke.restrictions),\n",
       " EntityKey(cause.diabetes_mellitus_type_2.restrictions),\n",
       " EntityKey(cause.asthma.restrictions),\n",
       " EntityKey(cause.gout.restrictions),\n",
       " EntityKey(cause.osteoarthritis.restrictions),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_hypertension.restrictions),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_glomerulonephritis.restrictions),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_other_and_unspecified_causes.restrictions),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_diabetes_mellitus_type_2.restrictions),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.exposure),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.distribution),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.exposure_standard_deviation),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.exposure_distribution_weights),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.population_attributable_fraction),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.tmred),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.relative_risk_scalar),\n",
       " EntityKey(population.age_bins),\n",
       " EntityKey(healthcare_entity.outpatient_visits.utilization),\n",
       " EntityKey(healthcare_entity.inpatient_visits.cost),\n",
       " EntityKey(healthcare_entity.outpatient_visits.cost),\n",
       " EntityKey(cause.ischemic_heart_disease.cause_specific_mortality),\n",
       " EntityKey(cause.ischemic_stroke.cause_specific_mortality),\n",
       " EntityKey(cause.diabetes_mellitus_type_2.cause_specific_mortality),\n",
       " EntityKey(cause.asthma.cause_specific_mortality),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_hypertension.cause_specific_mortality),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_glomerulonephritis.cause_specific_mortality),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_other_and_unspecified_causes.cause_specific_mortality),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_diabetes_mellitus_type_2.cause_specific_mortality),\n",
       " EntityKey(cause.ischemic_heart_disease.prevalence),\n",
       " EntityKey(cause.ischemic_heart_disease.disability_weight),\n",
       " EntityKey(cause.ischemic_heart_disease.excess_mortality),\n",
       " EntityKey(cause.ischemic_stroke.prevalence),\n",
       " EntityKey(cause.ischemic_stroke.disability_weight),\n",
       " EntityKey(cause.ischemic_stroke.excess_mortality),\n",
       " EntityKey(cause.diabetes_mellitus_type_2.prevalence),\n",
       " EntityKey(cause.diabetes_mellitus_type_2.disability_weight),\n",
       " EntityKey(cause.diabetes_mellitus_type_2.excess_mortality),\n",
       " EntityKey(cause.asthma.prevalence),\n",
       " EntityKey(cause.asthma.disability_weight),\n",
       " EntityKey(cause.asthma.excess_mortality),\n",
       " EntityKey(cause.gout.prevalence),\n",
       " EntityKey(cause.gout.disability_weight),\n",
       " EntityKey(cause.osteoarthritis.prevalence),\n",
       " EntityKey(cause.osteoarthritis.disability_weight),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_hypertension.prevalence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_hypertension.disability_weight),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_hypertension.excess_mortality),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_glomerulonephritis.prevalence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_glomerulonephritis.disability_weight),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_glomerulonephritis.excess_mortality),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_other_and_unspecified_causes.prevalence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_other_and_unspecified_causes.disability_weight),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_other_and_unspecified_causes.excess_mortality),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_diabetes_mellitus_type_2.prevalence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_diabetes_mellitus_type_2.disability_weight),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_diabetes_mellitus_type_2.excess_mortality),\n",
       " EntityKey(cause.ischemic_heart_disease.incidence),\n",
       " EntityKey(cause.ischemic_stroke.incidence),\n",
       " EntityKey(cause.diabetes_mellitus_type_2.incidence),\n",
       " EntityKey(cause.asthma.incidence),\n",
       " EntityKey(cause.gout.incidence),\n",
       " EntityKey(cause.osteoarthritis.incidence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_hypertension.incidence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_glomerulonephritis.incidence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_other_and_unspecified_causes.incidence),\n",
       " EntityKey(cause.chronic_kidney_disease_due_to_diabetes_mellitus_type_2.incidence),\n",
       " EntityKey(risk_factor.high_body_mass_index_in_adults.relative_risk)]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "art.keys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def format_for_dismod(df, data_type):\n",
    "    df = df.query('draw==0 and sex==\"Female\" and year_start==2017').copy()\n",
    "    df['data_type'] = data_type\n",
    "    df['area'] = 'all'\n",
    "    df['standard_error'] = 0.001\n",
    "    df['upper_ci'] = np.nan\n",
    "    df['lower_ci'] = np.nan\n",
    "    df['effective_sample_size'] = 10_000\n",
    "\n",
    "    df['sex'] = 'total'\n",
    "    df = df.rename({'age_group_start': 'age_start',\n",
    "                                     'age_group_end': 'age_end',}, axis=1)\n",
    "    return df\n",
    "\n",
    "p = format_for_dismod(art.load('cause.ischemic_heart_disease.prevalence'), 'p')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "i = format_for_dismod(art.load('cause.ischemic_heart_disease.incidence'), 'i')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = format_for_dismod(art.load('cause.ischemic_heart_disease.excess_mortality'), 'f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "m_all = format_for_dismod(art.load('cause.all_causes.cause_specific_mortality'), 'm_all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "csmr = format_for_dismod(art.load('cause.ischemic_heart_disease.cause_specific_mortality'), 'csmr') # could also try 'pf'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "dm = dismod_mr.data.ModelData()\n",
    "dm.input_data = pd.concat([p, i, f, m_all,\n",
    "                          csmr\n",
    "                          ], ignore_index=True)\n",
    "for rate_type in 'ifr':\n",
    "    dm.set_knots(rate_type, [0,40,60,80,90,100])\n",
    "\n",
    "dm.set_level_value('i', age_before=30, age_after=101, value=0)\n",
    "dm.set_increasing('i', age_start=50, age_end=100)\n",
    "\n",
    "dm.set_level_value('p', value=0, age_before=30, age_after=101)\n",
    "\n",
    "dm.set_level_value('r', value=0, age_before=101, age_after=101)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "m_all    23\n",
       "f        23\n",
       "p        23\n",
       "csmr     23\n",
       "i        23\n",
       "Name: data_type, dtype: int64"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dm.input_data.data_type.value_counts()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/ihme/homes/abie/.local/lib/python3.6/site-packages/pandas/core/indexing.py:480: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[item] = s\n",
      "/ihme/homes/abie/.local/lib/python3.6/site-packages/pandas/core/indexing.py:480: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame.\n",
      "Try using .loc[row_indexer,col_indexer] = value instead\n",
      "\n",
      "See the caveats in the documentation: http://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n",
      "  self.obj[item] = s\n"
     ]
    }
   ],
   "source": [
    "dm.setup_model(rate_model='normal', include_covariates=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pymc as pm\n",
    "m = pm.MAP(dm.vars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Current log-probability : 1347.572455\n",
      "Current log-probability : 1377.997951\n",
      "Current log-probability : 1396.019187\n",
      "Current log-probability : 1400.324203\n",
      "Current log-probability : 1400.952332\n",
      "Current log-probability : 1401.056619\n",
      "Optimization terminated successfully.\n",
      "         Current function value: -1401.056619\n",
      "         Iterations: 6\n",
      "         Function evaluations: 7548\n",
      "CPU times: user 26.7 s, sys: 43.6 ms, total: 26.7 s\n",
      "Wall time: 26.7 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "m.fit(verbose=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.core.pylabtools import figsize\n",
    "figsize(11, 5.5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAqAAAAFaCAYAAADB3e8nAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeZgU1dn38e89OzMgIKCPDoxgQONGREbcYp4kasCVxGBETUTEkBiNJiYxmkXFJWpi1LiRENe4K3HhcQluGPOqKCAKglERURYXZIfZp+/3j65peoaevbure+b3ua6+pk/VqaqbhinuPufUOebuiIiIiIikS07YAYiIiIhI96IEVERERETSSgmoiIiIiKSVElARERERSSsloCIiIiKSVkpARURERCStlICKiIhI0pnZIjP7ethxSGYyzQMqIiIiIumkFlARERERSSsloJLVzGyZmV1oZovNbJ2Z3WFmRWHHJSLS3QX358PDjkMykxJQ6QpOAUYDXwJ2A34XbjgiIiLSEiWg0hXc5O7L3X0tcAVwUtgBiYiISPOUgEpXsDzu/UfAzmEFIiIiIq1TAipdwaC492XAqrACERERkdYpAZWu4CwzG2hm2wO/AR4MOyARERFpnhJQ6QruA54Blgavy8MNR0RERFqiieglq5nZMuAMd38u7FhERESkbdQCKiIiIiJppQRURERERNJKXfAiIiIiklZqARURERGRtFICKiIiIiJplRd2AMnUv39/Hzx4cNhhiEgI5s2b94W7Dwg7DmlM92WR7FO7YEGHjssfPrxRuaX7ckoTUDMbA/wFyAVudfermuz/GnA9MBwY7+7T4/bVAwuD4sfuflxr1xs8eDBz585NVvgikkXM7KOwY5Bt6b4skn1Wlg5qvVICpU1+11u6L6csATWzXOBm4AhgBTDHzGa4++K4ah8DpwG/THCKSnffN1XxiYiIiEg4UjkGdBSwxN2XunsN8AAwNr6Cuy9z9wVAJIVxiEgKRTZupO7jj/GIfo1FRKRtUpmAlgLL48orgm1tVWRmc81stpl9O7mhiUhn1b7/PuvO/zWfjNiPzw46hM8O/iobr7ue+rVrww5NREQyXCrHgFqCbe2ZdLTM3VeZ2a7AC2a20N0/2OYiZpOByQBlZWUdi1RE2qz+k0/Y+Mc/UfHwdIibR7h++XI2XfNnttxzLwMef5S8gQNDjFJERDJZKltAVwDxo1gHAqvaerC7rwp+LgVeBEY0U2+au5e7e/mAAXoAViRVvL6eTX/9K58d+r9UPPRwo+STwsLY28inn7LmBxOIrF8fQpQiIpINUpmAzgGGmdkQMysAxgMz2nKgmfU1s8LgfX/gEGBxy0eJSKrULVvGF+NOYONlV+CVlbHthd/8Jv0fmc7Oi9+mz7XXQH5+tP5777Fm0hl4dXVYIYuISAZLWQLq7nXA2cBM4B3gIXdfZGaXmtlxAGa2v5mtAE4A/mZmi4LD9wDmmtlbwCzgqiZPz4tImlQ++xyfjzmKmtfnxLbl7fFl+t1/H/3vvovCAw7AioooOfFE+l7351idmtmvsWnqX8MIWUREMlxK5wF196eAp5psuyju/RyiXfNNj3sF2CeVsYlIy9ydzTfexMY//mlrd3tuLr1+di69fno2FrR2xiv+zneoX76CjVf/EYDNN91MyYnfI3enndIZuoiIZDgtxSkiVNfWNypXVdWw/vxfRxPJIPnMLS1lwIzH2O68n1PT5NYRf3zPs35C3h57AOCVlWz4Q6P1J0RERLrWUpwi0jGF+bkcePFMAPLq65jxxRNU/t8Tsf0fDd6Lf570Cyrm18D815k6cVSsPsDsKaNj7y03lz5TLuGL750IQOUjj1Bz2gQKRu6Xpj+NiIhkOrWAikhMbqSOXzz/10bJ54J9v869Ey+ioqR3m89TeMjBFB11ZKy84Uq1goqIyFZqARURAHIiEX7679sZ9fGbsW0lp09kzJRLODKn/d9Ve//ut1Q98yzU1VHz6qvUzJ9PwYiEs6mJiEiG8KoqSlcub71iM8daUVGb6qoFVESoqqnjyZqXOPSD12Pbev7kTHpfOgVLkHxW19Yze8ro2KvpGFKAvF12ocfYravvbrpFT8SLiGS6hgRyyz33srJ0ECtLB7H628e369i2UAIqItTe+ncq7rknVi459Qds95sLMUu0oFl0zGhL5Qa9zvxR7H3V009Tt/TDJEQrnWVmY8zsXTNbYmYXJNhfaGYPBvtfM7PBcfuGm9mrZrbIzBaaWdv/xxGRrFH5xJOx9z2OOTrp51cCKtLNVT79NBv/cGWs3OP479D7isubTT7bI3+PPSj85jeiBXc2/W1ap88pnWNmucDNwJHAnsBJZrZnk2qTgHXuPhS4Drg6ODYPuAf4sbvvBXwdqE1T6CKSJvVr11L9yiuxco+jjkr6NZSAinRRTbvFE3WT17y9iHU/PTc21VLBgQfQ98/XUFPv29TtqF5nnhl7X/Hww9SvXZu0c0uHjAKWuPtSd68BHgDGNqkzFrgreD8dOMyi30i+BSxw97cA3H2Nu2/7D0tEslrV0/+C+uivdkF5Obk7J38uZz2EJNJFxU+tBNGpks68Y+sYz8LKLZxxy6/oGyyt+cl2A7hw6Hg2XTGr0bRKnVVw0IHkD9+H2gULobqaigcepNdPzmz9QEmVUiD+CYMVwAHN1XH3OjPbAPQDdgPczGYCA4AH3P2PTS9gZpOByQBlZWVJ/wOISGpVPrF1JpRUdL+DWkBFuid3jnvkRvqu+wwA69mTK791DpuKeiX9UmZGyWmnxcpb7r4Hr1ejWYgSja1o2uTdXJ084KvAKcHP75jZYdtUdJ/m7uXuXj5gwIDOxisiaVS/di3VL2/tfi9KQfc7KAEV6VamThzF1Imj+GPNfHZ/Z+va7n3/fA0r+6Ruuczi447F+vQBoP7jj6me9WLKriWtWgEMiisPBFY1VycY99kbWBts/7e7f+HuFUSXWtYKAyJdSNXMZ2Ld7/n77Ude6c4puY664EW6qIapkuLLhfm51CxcyMartvaalpwxiZzRY5h9TO42dZPFevSgZPyJbP7r3wDYfNc/KDp8m4YzSY85wDAzGwKsBMYDJzepMwOYALwKjANecPeGrvfzzawYqAH+l+hDSiLSRaSj+x3UAirSZSWaKilSWcm6s8+B2uiDy/kj9qX3b3/T5mmVOqPkB9+H4Mn66lmzqPvoo6RfQ1rn7nXA2cBM4B3gIXdfZGaXmtlxQbXbgH5mtgQ4D7ggOHYdcC3RJPZN4A13f7LpNUQkO0XWraP6/70cK6cyAVULqEg3svGyy6lbsgQAKy5m+xtvwAoK0nLtvMGDKfzG16l+YRa4s+Xe++j9mwvTcm1pzN2fItp9Hr/torj3VcAJzRx7D9GpmESki6mcORPq6gDIHzGCvNLSlF1LLaAi3UTVS/9hy13/iJV7XzqFvCFD0nLthimgSn7wg9i2igcexGtq0nJ9ERFpXaonn4+nFlCRbiCyZQvrz/91rFw0ZjTF409M2/UbpoTKidQztaQv/besI7JmDTddcDOL9zkkbXGIiEhikfXrqf7P/4uVU52AqgVUpBvYeNXV1C+PTv1ofXrT58o/JGWlo/aK5OTy3O6Hxsr7vT6zhdoiIpIujbrf9/0KeQMHpvR6agEV6eKq58xhy+13xMp9pkwhd4cdQovn+d0PZfxbT0J9PYM/XMQNh25P/tChnT7vX09PQnAiIt1U4+73Y1J+PSWgIl2Y19ay/oKtD/oUfvOb9Pju8WmPo+mUUKu/eIaaZ58FohPT95lySdpjEhGRqG26349OzeTz8dQFL9KFbb71Nur++y4Qfeq9z1XJ6XpvyzrzLSk65ZTY+4qHpxMJlgMVEZH0q3zm2a3T831lOHlpWEJXLaAiXVTdihVs+vO1sXKvX/w8aVNqtLbOfFNTJ45qVP/Vi49gy+BdqF/2Eb5hA5WPPkbJySclJTYREWmfdHe/g1pARbqsDb+/CA9aFvP2+DI9J00KOaKtLCen0ZRMW+68C/emy5GLiEiqRTZsoPqll2LlVD/93kAtoCJdUOXMmVQ982ys3OfKP2D5+Sm95tSJo9pVv+TE77HxT3+CqmpqFy2iZu48CvcvT1F0IiKSSFV89/vwfdLS/Q5KQEW6nMiWLWz4XWxRG4pPOZnC/fdP6jWaW2e+XfX79qX429+m4oEHAdhy111KQEVE0iyM7ndQF7xIl7Ppz9dSv2oVADn9+tH7wguSfo32rh3fXP2S0ybEtlU+8ST1n3ySpAhFRKQ1kY0bqYrvfk/D0+8NlICKZIm2PHle8/YiNt96W6zc+6Lfk9O3b7P1w1awzz4UNLTO1tay+Y47Q41HRKQ7qXr2OQiWRM7fZx/yBg9O27XVBS+SJVp78twi9Zw27beU1kcTzYKDD+abC0pg4cxY/UzU88eTWTtnDhCdE7TXOT8lp2fPkKMSEen6Kp94IvY+na2foBZQkS5jv9efoXTF+wDU5ebR56orIYTlNtur6FvfIm/XXQHwjRvZcu99IUckItL1RTZtourf6X/6vUHKE1AzG2Nm75rZEjPbZjCamX3NzN4wszozG9dk3wQzez94TWh6rEh3N3XiKKZOHMVNYwZx1IsPxLb3/fm55H9p1xAjazvLyaHnjybHyltuvQ0PnsgUEZHUqHr2OaiuBiB/r73IGzIkrddPaRe8meUCNwNHACuAOWY2w90Xx1X7GDgN+GWTY7cHLgbKAQfmBceuS2XMIpmqpSfP1190Mb55MwB5X/oSvX5yZrufVA9T8bjvsvFP1xD54gvqV62i4pFHKDnxxLDDEhHpsiqfjH/6Pb2tn5D6FtBRwBJ3X+ruNcADwNj4Cu6+zN0XAJEmx44GnnX3tUHS+SwwJsXximSs5p4kr5w5k6qnno5t73PVlVhhYbufVA+TFRXRc9LpsfKma6/Hg2/mIiKSXJHNm6ma9WKsnM7plxqkOgEtBZbHlVcE25J2rJlNNrO5ZjZ39erVHQ5UJBtFNm1i/W9+FysXnzSewoMPCjGijis5fSI5228PQP2KFWy5//6QIxIR6Zqqnovrft9zT/J2TW/3O6Q+AU30BERb19tr07HuPs3dy929fMCAAe0KTiTbbbz6j0Q+/RSAnP796f3b34QcUcfl9OxJz7PPipU3/eVGIsFSoiIikjyNJ59Pf/c7pD4BXQEMiisPBFal4ViRLq967jy23HlXrNx7ysWxOT+zVc9Tf0DO/+wIQOTzzxv9+UREpPMiW7ZQNWtWrFwUQvc7pD4BnQMMM7MhZlYAjAdmtPHYmcC3zKyvmfUFvhVsE+n2vKaG9eefDx7tFCj85jfoMXZsK0dlPuvRg+3OPTdW3nTTzUQ2bQoxIhGRrqXqueegKtr9nrfHHqHNmJLSBNTd64CziSaO7wAPufsiM7vUzI4DMLP9zWwFcALwNzNbFBy7FriMaBI7B7g02CbS7W26+Rbq3n0PACsups+Vf8CyYM7PtigefyK5ZWXkDhlCnysuw0pKwg5JRKTLyITud0jDSkju/hTwVJNtF8W9n0O0ez3RsbcDt6c0QJEsU/vee2y64cZYebsLfk3ewIS/QlnJCgrof9895A4ahOVpsTYRkWSJbNlC1QsvxMphJqBaCUkki3h9Pet/ef7WtXtH7EvJaV1vjYa8IUOUfKZQGxYIKTSzB4P9r5nZ4Cb7y8xss5n9sumxIpK5qp57Ptb9nrv77uQPHRrbV11bn9ZYdIcXySJb7ryLmnnzooX8fPpe8ycsN3Pn95TM08YFQiYB69x9qJmNB64G4lcGuA54GhHJKvHd78XHHsOBF299tGb2lNGcecfraYtFLaAiWaLu44/ZeOVVsXKvc35K/pe/HGJEkqVaXSAkKDdMQTAdOMyCQcZm9m1gKbAoTfGKSBJEKiqozpDud1ALqEhWcHfW//J8PJgXM+/Lu9Mrbs5MkXZItMjHAc3Vcfc6M9sA9DOzSuDXRFtPm+1+N7PJwGSAsrKy5EUuIh1W9dzzeFUVAHm770b+sGFEv0tuNXXiqKRe86+nN79PCahIFqi49z6qX345WsjJiXa9FxSEG5Rkq7Ys8tFcnSnAde6+uaVZF9x9GjANoLy8vK2Lj4hICjV++v0YqmvrmT1ldGxbdW19WpdsVgIqkuHqVq5kw2WXx8o9fzSZghEjQoxIslxbFvloqLPCzPKA3sBaoi2l48zsj0AfIGJmVe5+U+rDFpGOilRUUP3887Fyj2OOJr9JspnO5BOUgIpkNHdn/a/OxzdvBiBv113Z7hfnhRyVZLnYAiHASqILhJzcpM4MYALwKjAOeMHdHTi0oYKZXQJsVvIpkvmqX5i1tft92DDyd9st5IiUgIpktIr77qf63y9FC2b0ufYarEePcIOSrBaM6WxYICQXuL1hgRBgrrvPAG4D7jazJURbPseHF7GIdFblE0/E3of98FEDJaAiGapuxQo2XHpZrNzzh2dQuP/+IUYkXUUbFgipIro6XUvnuCQlwYlIUkUqK6PzfwYyJQHVNEwiIWk66W982SOR6FPvcV3vhef9osXjRUREmqp+YdbWGVSGDiVv991DjihKLaAiISnMz212EuCRs5/myP/8B4CI5fD3w8/gd71KtqkvIiLSkqbd7y3NYJFOagEVyTDbf7GKw2f+I1aefcixrCzLjG+sIiKSPTxDu99BLaAiGeWWH+zH6u9cQW1tdK33vC/vzrh//JkTCgtDjkxERLJN1awX8YoKAPK+9CXyMmj1PCWgIiFJNAlw9c03U/vGG9EN+fn0/ctfsCD5DHvSYBERyS6VT8ZNPn/0URnT/Q7qghcJTdPk0RYtZNN118fK2533cwr23qvZ+ko+RUSkOV5ZSdWzz8XKPY45JsRotqUEVCQDRCoqWPfTc6GuDoCC/fen51k/CTkqERHJVlX//je+ZQsAuUOGkLfnHiFH1JgSUJEMsGHKZdQtXQqAlZTQ94brsVy1cIqISMc0Xvs9c55+b6AEVCRklf/6FxX33BMr977sUvLKykKMSEREsplXVWV09zsoARUJVf0nn7DuF7+KlYuOPpri77W4AI2IiMg24hcnqXrppdhCJrmDB5O/154t1k9UTjU9BS8SEq+vZ+05P8PXrwcgd+ed6fvHqzKum0RERDJf/OIm57x4K/8bbH+pbD9m3Tlnm/pTJ44KdXETtYCKhGTTDTdS88or0UJODn1vuoGcPn3CDUpERLJaXn0t5R+9FSu/s/dBIUbTPLWAioSgevZsNl17Xazc69xzKDzggBAjEhGRbOFVVVhR0TbbY62Yl28d83lxG49NNyWgImlWv3Yta8/6KUQiABQceAC9fnZuyFGJiEi2sKIiVpYO6tCxpSuXA+EvbqIueJE08vp61p39UyKffgpATt++bH/TjVievguKiEj6hL24iRJQkSRq7anCTTfcSPW/X4qV+/7lenJ32ilhXRERka5KzS4iSRT/FCJEx+OcecfrAAxZ8hYn33UtDc+49/zp2Xz9pVp4aWasroiISHegFlCRNOi97nO+89B1mDsAy4bsxXa//EXIUYmIiIQj5QmomY0xs3fNbImZXZBgf6GZPRjsf83MBgfbB5tZpZm9Gbz+mupYRVLhlvH7cN4zt1BcsQmAnB134MBH/qFxnyIi0qywJ4pPtZT+D2hmucDNwBHACmCOmc1w98Vx1SYB69x9qJmNB64GTgz2feDu+6YyRpFkavpUYVVNHZW/+S21CxdGN+Tns/3f/kbuDjuE/gSiiIhkrpaGdE2dOCqssJIm1S2go4Al7r7U3WuAB4CxTeqMBe4K3k8HDjMtBSNZqmkCWXfXnVQ89HCs3PuSiyncvzxhXSWfIiLSXaQ6AS0FlseVVwTbEtZx9zpgA9Av2DfEzOab2b/N7NAUxyqSVFUvvsiGSy+LlYu/dwIlE04NMSIREQlTZ7rV3Z3rvlTNZU9fm+ywQpHqQWiJWjK9jXU+AcrcfY2ZjQQeM7O93H1jo4PNJgOTAcrKypIQskjn1S5Zwtozz4pNNp+/3370ufIPWuddMoKZjQH+AuQCt7r7VU32FwL/AEYCa4AT3X2ZmR0BXAUUADXAr9z9hbQGL5LFWupWb+r6749k9pTReCRC1b9m8vlRx1C3YEG6Qk25VCegK4D4qfoHAquaqbPCzPKA3sBad3egGsDd55nZB8BuwNz4g919GjANoLy8vGlyK5J29WvWsOYHE/CN0e9KuTvtRL/b/p4RS5+JdHJs/hfAse6+ysz2Bmayba+WiNCG5TIDicZzbli3iQKv58KTpvDtt55m4IZPUxZnWFKdgM4BhpnZEGAlMB44uUmdGcAE4FVgHPCCu7uZDSCaiNab2a7AMGBpiuMV6RSvqmLt6WdQ//HHAFiPHmx/+63k7rBDyJGJxMTG5gOYWcPY/PgEdCxwSfB+OnCTmZm7z4+rswgoMrNCd69Ofdgi2aWzy2V+uv8BnL2qSZtdYSEl409MfFCWSWkC6u51ZnY20W/JucDt7r7IzC4F5rr7DOA24G4zWwKsJZqkAnwNuNTM6oB64MfuvjaV8Yp0hkcirDv359TMDRrpzeh7840UDB8ebmAijSUam39Ac3WC+3jD2Pwv4up8F5iv5FMkNerjkk/r1Yui73+f3pPP6DINGimfiNDdnwKearLtorj3VcAJCY77J/DPVMcnkgzuzoaLL6HyiSdi23pf9Ht6jNbqRpJxOjM2P7rTbC+i3fLfSngBjc0XSYqc/v3pecYkSiacSs5224UdTlJpJmyRJNh88y1suf2OWLnk9ImU/PCMECMSaVaHx+YDmNlA4FHgVHf/INEFNDZfpPN6X3EZJSeeiPXosc0+r6qidOXyBEe1rrmxqemmpThFOmnLP+5m45VbHyLuccwx9J5yiZ54l0wVG5tvZgVEhz3NaFKnYWw+NB6b3wd4ErjQ3V9OW8Qi3VDP005LmHwCnUogMyH5BCWgIp1S8cijrP/Nb2PlgoMOou9frsNy9KslmSmYb7lhbP47wEMNY/PN7Lig2m1Av2Bs/nlAwzLKZwNDgd/HLZPcNQakiXRAV18uM5XUBS/SgqbLY8aXK594knU/+zl4tIcxf9+v0HPatEbfLrW8pmSiTozNvxy4POUBimSJrr5cZiopARVpQXM3ly+//SrHP3QtOcFE85/vWMbdR/2ca7fvs019ERHpPvp/9jHD578ISkBbpARUpJ32XPgy3374+ljy+cWAUu497SIqi3uFHJmIiIShZ9VmDv3gdT4/6gZ+/FbDakU3hRpTplMCKtIOW+67n+Mfui7W7Z43dCh7P/QA1+y4Y8iRiYhIOnltLZuefZ4nVz1C1bPPQW0ttWEHlUWUgIq0oLq2ProWrzub//Y31v/qiti+vGHD6P/g/eTGJZ8N9ePLGgMqItJ11C5azJaHHqLy0ceIrFmzbYWCAnp8K+EUuRJHCahICwrzc/H6ejZcMqXRPJ/5w/eh3733kLv99tvUb6ksIiLZp371aioffYyKh6dTu3hxwjr5I0ZQfMI4io87lpy+fdMcYfZRAirSgkhFBet+eg5V/9r6YFHBAaPod+cdXW5VChER2cqrq6l67nkqHnqYqlmzoH7bKZZy/mdHiseNo/iEceQPHRpClNlLCahIM+pXfcKa0ydRu3BhbFuPY4+h7/XXZcxEviIikjzuTu1bb1Hx8HQqHnscX79+20pFhfQ48kiKTxhH4Ve/iuWqp6sjlICKJFAzfz5rJp1B5LPPY9t6/mgy2/3ut5pkXkSki6n/9FMqHnmUioenU/feewnrFIzan+ITTqDHMUe3qQesKyyXmUpKQKXbaWlyeYAtDzzA+gt/CzU10Q15efS54nJKvn9KtK7yTxGRrOeVlVQ+8wwVD0+n+t8vQTC1Xrzc0lKKx32X4hPGkTdkSLvO3xWWy0wlJaDS7TQ3uXxubQ3feuoORs55JrbP+vTmooN/yNvv94eLZ2pieRGRLObu1MydR8XD06n8v//DN27cpo4VF1N01FGUfO8ECg46UL1eKaIEVATou+ZTjn/gGnb65MPYts92LGP4o/fx9u3/DTEyEZHuqzNd0fHH1q1cSeX0f7Ll4enUf/hhwvoFBx1E8Qnjol3sJSUdjlnaRgmodHsVjz3G2bf+Bt+0Kbatx3HHsu+fryGnuBhQAioiEgYrKmJl6aAOHVu6cjkV/3yEiocepvrll2MLiMTLHbxL9Cn2cd8lb1D0OtW19RTG1WlpPufWhnRJ85SASrfTMFl8ZNMm1v/uItadNX3rzvx8el98ESWnTcDMNLG8iEgWW3fOudtss1696HHsMRSfMI6C/ffHzBrtb26YViJTJ47apq60jRJQ6XYK83Op/n8vs+68X1C/cmVse+4uZWx/y80U7Ltvo7pNjxURkSxjRuHXDqX4hHEUjRlDTo8eYUfU7SkBlazXni6QyKZNrLviSqruvrvR9h7f/S59rriMnF69UhqriIikT97QodHViY4/ntydd+rweaZOHJXEqASUgEoX0KbuEnd2e+d1xjxxK9ttXBvbvKmwhLLrrqZ47Nh0hSsiImmyw4svbNPF3pr2DL3SMK2OUwIqXV7fL1Yx+qnbGfre/Ebb5w4aztRDT+VfSj5FRDKKu1P/0UfkDR7cqfO0N/mE9g290jCtjlMCKl3S1ImjiKxfz6YbbmTzHXdunVQeyOnfnz/tczyv7Lo/dODmJCIiyRXZsoXatxZQM29e9PXGfCJr1nR4JSHJfEpAJeO0d1qLpl0gVRs3U3vvPWy66SZ8/YatFc0oPuUUin75S64d0K/N5xcRkeRxd+o/XBYkmm9QM+8Nat95J+FKRMnQkf9TNLVS6ikBlYzTnikwYOs0GEW1VRz+35eYtHQWkdWrG9UpGDmS3pdfSsHw4QmvJyIiqRHZvJma+W9S+8YbVM97g9o33iCybl2rx1nv3km5fkf/T4mvL8mnBFSyXv1nn3HivMc5cvEL9KreQvx36Nxdyuj9m99QdPRRHRoLJCLJl6zVbaTt0vWZeyRC3dKl1MyLtmzWvDGPuv++m3AS+EbMyNt9NwpGjqRgvxHwlREU7z6sQ/HGO/OO1/UEe4ZSAipp0dkujaY3EI9EqJn9GlvuvZdPp1wOhCQAACAASURBVDzF92prG+3P3Wknep17DsUnfg8rKOhc8CKSVJ1d3UbaL1WfeWTjRmrefDNIOOdRM39+46FPzcXTpw8F++1HwX4joknniH05+JpXojvfgdnjd09p8qikNHxKQCUt2tMFcv33RzY7rUXtu+9S+X9PUPHPR6j/+ONtjs3dpYyiM35I75PHq5VERCSJPBKhbsmSrcnmG29Q9977rbdu5uSQu9tuFJaXxxLOvC/tmpZeqakTR7V7qiRNrZQeKU9AzWwM8BcgF7jV3a9qsr8Q+AcwElgDnOjuy4J9FwKTgHrgHHefiWSEVA7Sjk9WC+pqeGZMXzbMmkXV8y9Qt2RJwmMK9t+fkomn0ePoo7A8fa8SaYnuy9Jen48/mbq33sI3bmy1bk7fvuTvtx93rO/Fuzt+iSX9h/DilWOjjQ4VwH/WRF9xErVIhtVKqamV0iOl/1ObWS5wM3AEsAKYY2Yz3H1xXLVJwDp3H2pm44GrgRPNbE9gPLAXsDPwnJnt5u71qYy5u2pvQtmRQd2tbfPqauqWLaPiscc49bUnGPb5Uoat/pA1d9YlPKf17k3xt8dSfPJJFOy9d7PXFpGtdF/uulLZMFD7n/8k3B7JyeHzHXdhxaDdWDloN1aU7c667f+HqacfwD8vbvt3k1S2OiqhzEypbioaBSxx96UAZvYAMBaIv9GNBS4J3k8HbrJou/xY4AF3rwY+NLMlwfle7WxQ62+9HVv9OfURJzdnaxdA03KibYnqpKt+qmP5238+jL2fcOgQHl/wSVDatntl9PCdOWXO0lh5w5Xz+caCVc2ee+3K55hZlBOczlnz+1fIqaoksn4Dkc8/p/7TT6lftSo2DUdzU8Nbjx4UHX4YPY45hqLDD1M3u0j7ZeR9uT3WXnFlt7kv19dHyM3Nabbc1LT/bL0vn3rorjwe3Je/d/8NzR7THuuLevHeDl/i8FPGcNPnxXxS+iVqC9p2H26tRTPRk+deVdXhcb96YC2zpToBLQXi/+WsAA5oro6715nZBqBfsH12k2NLkxFUzfTp1C5cmIxTdSnfjXu/+S04pIW6m1+C49tRv/KljseVN3QohV89hKJvfpOCgw8ip0ePjp9MRDLyvtwelbfcku5LZo3m78udS0D73ngDJ/1nC5/16g9mfOfs0VzYQv1kjbvsTAKp5DOzpToBTfS1rmlzWnN12nIsZjYZmAxQVlbW3vgkk5iRO2gQeUOHUrD3XuTvsw8F+5eTO2BA2JGJdCW6L0u75R57HI8fn7rJ2dVN3v2kOgFdAcTP+zAQaNpP21BnhZnlAb2BtW08FnefBkwDKC8vb+VRvKiS007l5vtf4czDhjH1+fdj2888LPGcY22pk676qTp3XcTJi+sGalpuumRlXX2EvLhuoKblVuXnk9OzJ9arJ7k77EDODjuSV7qzvrGKpF5G3pfbY7tfn98t7ssAmHHLc1vr/+TwFu7jnb0vJ5ESSmmNeWvTJ3Tm5NEb13vAYcBKYA5wsrsviqtzFrCPu/84GOx+vLt/z8z2Au4jOr5oZ+B5YFhLg93Ly8t97ty5rcbV8M2tLQO2U72EV3vqZ9K5RTKNmc1z9/Kw48h0mXJf7syclN3lvpzsWDT3qqRbS/fllLaABmOHzgZmEp3u43Z3X2RmlwJz3X0GcBtwdzCYfS3RJywJ6j1EdGB8HXBWsp60bPiFbMs3tPZ+i0tl/Uw6t4hkp0y9L7dHd7kvpzoWkTCltAU03draAioiXY9aQDNTKlpApWP0mUu6tXRf7lIJqJmtBj5qxyH9gS9SFE4qZWvckL2xZ2vckL2xtzfuXdxdT6xlmET35ZEjR46sXbCgQ+fLHz6cefPmzUtGbEmW0b9n+swzTrbGnrT7cpdKQNvLzOZmY4tJtsYN2Rt7tsYN2Rt7tsYtnZOtf+/ZGjdkb+zZGjdkb+zJjDucx+NEREREpNtSAioiIiIiadXdE9BpYQfQQdkaN2Rv7NkaN2Rv7Nkat3ROtv69Z2vckL2xZ2vckL2xJy3ubj0GVERERETSr7u3gIqIiIhImikBFREREZG06pYJqJmNMbN3zWyJmV0QdjwtMbNBZjbLzN4xs0Vmdm6wfXsze9bM3g9+9g071kTMLNfM5pvZE0F5iJm9FsT9oJkVhB1jImbWx8ymm9l/g8/+oGz4zM3s58G/k7fN7H4zK8rUz9zMbjezz83s7bhtCT9ji7oh+J1dYGb7hRe5pILuy+mj+3J66b6cWLdLQM0sF7gZOBLYEzjJzPYMN6oW1QG/cPc9gAOBs4J4LwCed/dhRNdjztQb9rnAO3Hlq4HrgrjXAZNCiap1fwH+5e5fBr5C9M+Q0Z+5mZUC5wDl7r430WUWx5O5n/mdwJgm25r7jI8EhgWvycDUNMUoaaD7ctrpvpwmui+3wN271Qs4CJgZV74QuDDsuNoR/+PAEcC7wE7Btp2Ad8OOLUGsA4N/rN8EngCM6AoKeYn+LjLlBWwHfEjwkF7c9oz+zIFSYDmwPZAXfOajM/kzBwYDb7f2GQN/A05KVE+v7H/pvpzWWHVfTm/cui838+p2LaBs/cfQYEWwLeOZ2WBgBPAasKO7fwIQ/NwhvMiadT1wPhAJyv2A9e5eF5Qz9bPfFVgN3BF0U91qZiVk+Gfu7iuBa4CPgU+ADcA8suMzb9DcZ5y1v7fSJln796v7ctrovhyelNyXu2MCagm2ZfxcVGbWE/gn8DN33xh2PK0xs2OAz909fv3gbPns84D9gKnuPgLYQoZ16yQSjMsZCwwBdgZKiHaRNJWJn3lrsuXfjnRMVv796r6cVrovZ55O/dvpjgnoCmBQXHkgsCqkWNrEzPKJ3uTudfdHgs2fmdlOwf6dgM/Diq8ZhwDHmdky4AGi3T3XA33MLC+ok6mf/Qpghbu/FpSnE73xZfpnfjjwobuvdvda4BHgYLLjM2/Q3Gecdb+30i5Z9/er+3La6b4cnpTcl7tjAjoHGBY8gVZAdDDwjJBjapaZGXAb8I67Xxu3awYwIXg/gegYpIzh7he6+0B3H0z0M37B3U8BZgHjgmoZFzeAu38KLDez3YNNhwGLyfDPnGgXz4FmVhz8u2mIO+M/8zjNfcYzgFODpy4PBDY0dAlJl6D7chrovhwK3ZebE/Zg15AG2B4FvAd8APw27HhaifWrRJu0FwBvBq+jiI7beR54P/i5fdixtvBn+DrwRPB+V+B1YAnwMFAYdnzNxLwvMDf43B8D+mbDZw5MAf4LvA3cDRRm6mcO3E90TFQt0W/Sk5r7jIl29dwc/M4uJPpEaeh/Br2S+u9B9+X0/hl0X05f3LovJ3hpKU4RERERSavu2AUvIiIiIiFSAioiIiIiaaUEVERERETSSgmoiIiIiKSVElARERERSSsloCIiIiKSVkpARURERCStlICKiIiISFopARURERGRtFICKiIiIiJppQRURERERNJKCaiIiIiIpJUSUBERERFJKyWgIiIiIpJWSkBFREREJK2UgIqIiIhIWikBFREREZG0UgIqIiIiImmlBFRERERE0koJqIiIiIiklRJQEREREUkrJaAiIiIiklZKQEVEREQkrZSAioiIiEhaKQEVERERkbRSAioiIiIiaaUEVERERETSSgmoiIiIiKSVElARERFJOjPb3czmm9kmMzsn7Hgks+SFHYCIiIh0SecDL7r7iLADkcyjFlDpcsxMX6xERMK3C7Ao7CAkMykBlS7BzJaZ2a/NbAGwRUmoiEh4zOwF4BvATWa22cx2CzsmySzm7mHHINJpZrYMWA8cC3zh7pXhRiQi0r2Z2YvAPe5+a9ixSOZRK5F0JTe4+/KwgxAREZGWqQteuhIlnyIiIllACah0JRpPIiIikgWUgIqIiIhIWikBFREREZG00lPwIiIiIpJWagEVERERkbRSAioiIiIiaaUEVERERETSSgmoiIiIiKRVl1oJqX///j548OCwwxCREMybN+8Ldx8QdhzSmO7LItmndsGCDh2XP3x4o3JL9+UulYAOHjyYuXPnhh2GiITAzD4KOwbZlu7LItlnZemgDh1X2uR3vaX7clK64M1sjJm9a2ZLzOyCBPsLzezBYP9rZjY42D7YzCrN7M3g9de4Y0aa2cLgmBvMzJIRq4iIiIiEq9MJqJnlAjcDRwJ7AieZ2Z5Nqk0C1rn7UOA64Oq4fR+4+77B68dx26cCk4FhwWtMZ2MVERERkfAlowV0FLDE3Ze6ew3wADC2SZ2xwF3B++nAYS21aJrZTsB27v6qR2fK/wfw7STEKiIiIiIhS0YCWgosjyuvCLYlrOPudcAGoF+wb4iZzTezf5vZoXH1V7RyThERERHJQsl4CClRS2bT9T2bq/MJUObua8xsJPCYme3VxnNGT2w2mWhXPWVlZW0OWkRERETCkYwW0BVA/ONSA4FVzdUxszygN7DW3avdfQ2Au88DPgB2C+oPbOWcBMdNc/dydy8fMEAzsIiIiIhkumQkoHOAYWY2xMwKgPHAjCZ1ZgATgvfjgBfc3c1sQPAQE2a2K9GHjZa6+yfAJjM7MBgreirweBJiFREREZGQdToBDcZ0ng3MBN4BHnL3RWZ2qZkdF1S7DehnZkuA84CGqZq+Biwws7eIPpz0Y3dfG+w7E7gVWEK0ZfTpzsYqIiLNa8OUemVmNisYt7/AzI4KI04RyX5JmYje3Z8Cnmqy7aK491XACQmO+yfwz2bOORfYOxnxiYhIy+Km1DuC6DCoOWY2w90Xx1X7HdFGhqnBdHtPAYPTHqyIZD2tBS8iItC2KfUc2C5435tmxuaLiLSmSy3FKSIiHZZoSr0DmtS5BHjGzH4KlACHpyc0Eelq1AIqIiLQtunvTgLudPeBwFHA3Wa2zf8jZjbZzOaa2dzVq1enIFQRyXZqARUREWjblHqTCJZFdvdXzawI6A98Hl/J3acB0wDKy8sTzuEsIpnJq6ooXbm89YrNHGtFRW2qqxZQERGBtk2p9zFwGICZ7QEUAWriFOlC4hPIzw4/gpWlg1hZOoiqWbPadWxrlICKiEhbp9T7BfDDYOq8+4HT3F0tnCJdUP2aNdS9899oIS+PglGjknp+dcGLiAjQpin1FgOHpDsuEUm/mldejb0vGDGCnJKSpJ5fLaAiIiIi0kj1K6/E3hcefFDSz68EVEREREQaqX45LgE9JPkdH0pARURERCSm/tNPqfvgg2ihsJCCkfsl/RpKQEVERES6mera+mbL1fHjP0eObNfT7W2lh5BEREREupnC/FwOvHhmrDx7ymjOvON1AI5+9HFGBNuf6VHG/wu2J5NaQEVEREQkZvDShbH3Hw3ZOyXXUAuoiIiIiDB14ijqli/ns99FFzezHj248LcnYQUFHTrfX09vfp8SUBEREZFuprq2ntlTRjcqF+bnNpp+qeCAUR1OPlujLngRERGRbqYwPzdhufrlrQ8gFR58cMqurwRURERERHB3ql9+OVYuPEQJqIiIiIikUP2Hy4h8+ikA1qsX+Xun5gEkSFICamZjzOxdM1tiZhck2F9oZg8G+18zs8HB9iPMbJ6ZLQx+fjPumBeDc74ZvHZIRqwiIiIisq1GrZ8HHoDlpe5RoU6f2cxygZuBI4AVwBwzm+Hui+OqTQLWuftQMxsPXA2cCHwBHOvuq8xsb2AmUBp33CnuPrezMYqIiIhIyxqv/5667ndITgvoKGCJuy919xrgAWBskzpjgbuC99OBw8zM3H2+u68Kti8CisysMAkxiYiIiEgbuXujFZBSsf57vGQkoKXA8rjyChq3Yjaq4+51wAagX5M63wXmu3t13LY7gu7335uZJSFWEREREWmi7r33iHzxBQA5ffuSt8eXU3q9ZCSgiRJDb08dM9uLaLf8j+L2n+Lu+wCHBq8fJLy42WQzm2tmc1evXt2uwEVEREQEql+Om//zoAOxnNQ+p56Ms68ABsWVBwKrmqtjZnlAb2BtUB4IPAqc6u4fNBzg7iuDn5uA+4h29W/D3ae5e7m7lw8YMCAJfxwRke6ptQdKgzrfM7PFZrbIzO5Ld4wikhqNxn+mcPqlBslIQOcAw8xsiJkVAOOBGU3qzAAmBO/HAS+4u5tZH+BJ4EJ3jz16ZWZ5ZtY/eJ8PHAO8nYRYRUQkgbgHSo8E9gROMrM9m9QZBlwIHOLuewE/S3ugIpJ0HolQ/Wp6JqBv0OkENBjTeTbRJ9jfAR5y90VmdqmZHRdUuw3oZ2ZLgPOAhm/WZwNDgd83mW6pEJhpZguAN4GVwN87G6uIiDSrLQ+U/hC42d3XAbj752mOUURSoHbxYnz9BgByBgwgb9iwlF8zKRM8uftTwFNNtl0U974KOCHBcZcDlzdz2pHJiE1ERNok0QOlBzSpsxuAmb0M5AKXuPu/0hOeiKRK/PjPwoMPIh3PfaduhlEREckmbXmgNA8YBnyd6Hj//5jZ3u6+vtGJzCYDkwHKysqSH6mIJFWjBDTF0y810FKcIiICbX+g9HF3r3X3D4F3iSakjejhUJHs4XV11Lz2WqxcePBBabmuElAREYG2PVD6GPANgOBB0d2ApWmNUkSSqnbBQnzzZgByd9qJ3MGD03JdJaAiItLWB0pnAmvMbDEwC/iVu68JJ2IRSYb46ZcKDjkkLeM/QWNARUQk0IYHSp3oTCbnpTk0EUmRxuu/p6f7HdQCKiIiItIteU0NNa/PiZXTMQF9AyWgIiIiIt1QzZtv4pWVAOTuUkbewIFpu7YSUBEREZFuqPH8n+lr/QQloCIiIiLdUuP5P5WAioiIiEgKeWUlNfPmxcpqARURERGRlKqZ9wbU1ACQN3QouTvumNbrKwEVERER6WaqX3459j7d3e+gBFRERESk26l+5dXY+3R3v4MSUBEREZFuJbJlCzVvvhkrFxx0YNpjUAIqIiIi0o3UvP461NUBkLfHHuT265f2GLp8AupVVaEc251l62eerXF39vqtHZut5xYRkcTCnP+zQZdfC96KilhZOqhDx5auXJ7kaLqHbP3MszVuSG3s2XpuERFJrNH6718NJwHt8i2gIiIiIhIV2bCB2oVvRws5ORQecEAocSQlATWzMWb2rpktMbMLEuwvNLMHg/2vmdnguH0XBtvfNbPRbT2niIiIiLRP9WuvQSQCQP4+e5PTu3cocXQ6ATWzXOBm4EhgT+AkM9uzSbVJwDp3HwpcB1wdHLsnMB7YCxgD3GJmuW08p4iIiIi0QyaM/4TktICOApa4+1J3rwEeAMY2qTMWuCt4Px04zMws2P6Au1e7+4fAkuB8bTmniIiIiLRDmOu/x0vGQ0ilQPzTACuApgMKYnXcvc7MNgD9gu2zmxxbGrxv7Zxp0dEHJLqzzj4cEtZnnq1xQ2pjz+Rzi4hI29WvXUvdO+9EC3l5FIwaFVosyWgBtQTbvI112rt924ubTTazuWY2d/Xq1S0GKiIizWvr2HszG2dmbmbl6YxPRDqnJm71o4KvfIWckpLQYklGAroCiG/iGAisaq6OmeUBvYG1LRzblnMC4O7T3L3c3csHDBjQiT+GiEj31dax92bWCzgHeC29EYpIZzWafinE7ndIThf8HGCYmQ0BVhJ9qOjkJnVmABOAV4FxwAvu7mY2A7jPzK4FdgaGAa8TbQFt7ZxpoS7C9MvWzzxb44bUxp7Nn0s3Ext7D2BmDWPvFzepdxnwR+CX6Q1PRDorfvxnQYgPIEESWkDdvQ44G5gJvAM85O6LzOxSMzsuqHYb0M/MlgDnARcExy4CHiJ6g/sXcJa71zd3zs7GKiIizUo0nr80voKZjQAGufsTLZ1IQ6NEMk/9Z59Rt2RJtFBQQGH5yFDjScpKSO7+FPBUk20Xxb2vAk5o5tgrgCvack4REUmZFsfem1kO0Wn0TmvtRO4+DZgGUF5ennD8voikV/WrceM/R+6H9egRYjRaCUlERKJaG3vfC9gbeNHMlgEHAjP0IJJIdmg8/dIhIUYSpQRUREQgbjy/mRUQHXs/o2Gnu29w9/7uPtjdBxOdQu84d58bTrgi0h6NHkA6+KAQI4lSAioiIm0dzy8iWahu5Urql30EgBUVUTBiRMgRJWkMaCbzqqoOP4XrVVVYUVGSI+r6svUzz9a4G66fqtiz9dzSfq2N52+y/evpiElEOq/R0++j9scKCkKMJqrLt4B25j8o/efWMdn6mWdr3J29fmvHZuu5RUQkKlPWf4/X5RNQERERke7K3al5JbMeQAIloCIiIiJZr7q2PmG5ftky6ldFJ7Swnj3JH75Pi/XTpcuPARURERHp6grzcznw4pmx8uwpoznzjtcZMedZjg62vVe6O5fd/QYAUyeO2qZ+OqkFVERERKSLGrx0Yez9sl33DjGSxtQCKiIiItIF3XLa/nz6l/eIBOWTfzae0/bOjCRUCaiIiIhIlquurW/UjV5dW0/Ohx8QWb0aAOvTm/w992yxfmF+btriVRe8iIiISJZrmjwW5uc2Xv3ooIOwnJwW66eTElARERGRLijT1n+PpwRUREREpIvxSITqV16NlTNh/fd4SkBFREREupjaxe/g69cDkNO/P3m77RZyRI0pARURERHpYhqtfnTwQZhZiNFsSwmoiIiISBeTyeM/QQmoiIiISJfidXVUv/ZarFx48MEhRpOYElARERGRLqR24UJ80yYAcv7nf8gdMjjUeBLpVAJqZtub2bNm9n7ws28z9SYEdd43swnBtmIze9LM/mtmi8zsqrj6p5nZajN7M3id0Zk4RURERLqLRk+/H3JIxo3/hM63gF4APO/uw4Dng3IjZrY9cDFwADAKuDguUb3G3b8MjAAOMbMj4w590N33DV63djJOERFphZmNMbN3zWyJmSW6n59nZovNbIGZPW9mu4QRp4i0rPrll2PvCw/JrOmXGnQ2AR0L3BW8vwv4doI6o4Fn3X2tu68DngXGuHuFu88CcPca4A1gYCfjERGRDjCzXOBm4EhgT+AkM9uzSbX5QLm7DwemA39Mb5Qi0hqvqaHm9TmxciY+gASdT0B3dPdPAIKfOySoUwosjyuvCLbFmFkf4FiiragNvht8y55uZoOaC8DMJpvZXDObuzpY71RERNptFLDE3ZcGjQIPEG1kiHH3We5eERRno0YDkYxT8+abeGUlALllZeQNzMxf01YTUDN7zszeTvAa29qxDadIsM3jzp8H3A/c4O5Lg83/BwwOvmU/x9ZW1m1P5D7N3cvdvXzAgAFtDElERJpotbGgiUnA0ymNSETarfH0S5n39HuDvNYquPvhze0zs8/MbCd3/8TMdgI+T1BtBfD1uPJA4MW48jTgfXe/Pu6aa+L2/x24urU4RUSkU1psLGhU0ez7QDnwv83snwxMBigrK0tWfCLSBo0S0AycfqlBqwloK2YAE4Crgp+PJ6gzE/hD3INH3wIuBDCzy4HeQKOn3BuS2qB4HPBOJ+MUEZGWrQDihzsNBFY1rWRmhwO/Bf7X3asTncjdpxFtXKC8vDxhEisiHedVVVhRUcJ9A6Y/1OFj06mzCehVwENmNgn4GDgBwMzKgR+7+xnuvtbMLgMaRsReGmwbSPQm9l/gjWCKgJuCJ97PMbPjgDpgLXBaJ+MUEZGWzQGGmdkQYCUwHjg5voKZjQD+RvRB0kQ9XiKSBlZUxMrSZh+PaVHpyuWtV0qDTiWgQVf5YQm2zyWuVdPdbwdub1JnBYm7fHD3CwlaSUVEJPXcvc7Mzibaa5UL3O7ui8zsUmCuu88A/gT0BB4OGg0+dvfjQgtaRLJWZ1tARUSki3D3p4Cnmmy7KO59s88EiIi0h5biFBEREZG0UgIqIiIikmGqa+tbLGc7dcGLiIiIZJjC/FwOvHhmrDx7ymjOvON1AKZOHBVWWEmjFlARERERSSsloCIiIiJZ4JYJ5fx5+0/DDiMp1AUvIiIikmGqa+uZPWU0AB6JsPHRx6m84Qbq3nuP4rHZP/uZElARERGRDFOYn4tHIlQ9+RQbr7uOunffCzukpFICKiIiIpJknVny0quqqHr+hWji+c5/G+2zkpJkhBc6JaAiIiIiSdbZ5TLXTv5R4/OVlFAy8TR6/uhHzRyVXZSAioiIdHGdbY3r6LHSeVZcTMnpEyk4fRI9dhwQdjhJowRURESki+tsa5yknxUX88jQr/H48NFsqurF7C6UfIISUBEREclQ3bnldsfZr3DPDXMT7vOqqg5/MciUz0UJqIiIiGSk7txym9uvX7P7OpNAZkLyCUpARUQkzbpzq5ZkL3fHKyuJrFuPr19PZP16qtesJXfTRiIbNhBZv57ateuw4P2Ahx/s1PXi5wFtKBfm53b2j5ExlICKiEhadedWrfZqmnS0lIS0p2535pEIvnEjkSCJbEgeI+vWx5LInI3BtvUbqF+3Dt+wgciGDVBTk7Y4m/7ddbW/SyWgIiIicVKZ9LW1vrtDTQ35FRUcdcW/KKqtobCuhjsnfIU/PTaf/Noa8muqKKipJq+2mvzaao7arS+V1VV4ZSVeUUH9lgqsKloe8PijHfkoYj772tehIB/Lz8fyC2LvycvHGt7nF8Te1+fmkVdYAPn5WEFBtFxUgOUXBHXzt/4syKcuN5/8ooLYeepycinoUUj+Hnt0Ku61Z51N3dp1EJdQ+saNEIl06rzSeUpARURE4hTm53LgxTNj5dlTRnPmHa83ruROXm0N14/bk2OvmElhXTWFdTXc/oPhXDPjLfJroklhfk0VBbXVsfIRQ/vw4MvvU1hXQ1FtNYeU9eLNjz4nv7aagoZjgldOkCT9Pe6yXzwGJzcTd+Vzyf0c4tV98EHqTt6CzrZ4Vz72eJIiiarLy6eiRy+qintSVVTCbnuU8eTSzWwpKGZTUU/O+u7+WJ8+5PTpk9TrdkWdSkDNbHvgQWAwsOz/t3fvMVKVZxzHv7/duezMogU06hawakpaaYNFNlWxsUZFwZhqU60YLyTFUBuN1pi2UO1FWxNMmgq2mlDZSAAACglJREFU3ohi1bTaFm2l1JQoimnSxAKxIl5Q6nUVRcVaUXaYHZ7+cd5lZ7czuyyzM3MO+3ySkznnnTNnn3ln5smz73nnDPBNM/ugwn5zgWvC5s/N7O7QvgboAHaE+041s62SssA9wHTgfeBcM3u1llidc84NTtIsYAnQCtxhZosG3J+Y3GylUjQSGEYDo9tofVdv+45PorZ++0S3V617hWzPTrLFAlvX3sS3t2yLishQJGaKBQDevg5uK/u77z0E5w0S18d/h5PLtgtvwMR6dIAbVHc2T3eunR35/diRG8OO3Bi6c2PozrUz+6tTWLSmi4+yY9iezXP7lTNpGTs2KixzuX7HKRRLXOTTHvZKrSOgC4DVZrZI0oKw/YPyHUKR+hOgEzBgvaQVZYXq+WY28DoD84APzOyzkuYANwDn1hirc865KiS1AjcDM4EuYG3I1c+V7RaL3PzeZZej7r5isfTJDrSjf3FJoVDT35hRtl58G5pyBcZ0GuVyKJ9DuTwtuRyWy9Gaz0dt+Txqa4tuczla8nmWPPEahVSWQirDtRccEx6bG/pvDeGgNY9xwZInSO3qIbWrxNILp3HTXzfSWirSUirRWurpt8zpnMCvVj5DaleJ1K4S35oxCYpFbOdOrNgDxZ1YTw/sLGI7C1hPD2tf2BLtXyrxhUPasWLt8y3HLVlMy9ixzF++ie3ZdrZn21l1/ZnRqf4qCsUS11/Sv6hMe1E54motQM8ETgzrdwNrGFCAAqcBj5jZNgBJjwCzgPuGOO5Pw/py4NeSZGZWY7zOOecq+zKw2cxeBpB0P1EuLi9AY5GbC3+qbT7jiGnLorao8IsKxTzW1kZre9jO5UNx2LZ7n1K2jfSY9t0FZE86S2a/MbTsLjRzuwtKpdPDno96zRX1GY1LT57Mawe8vHs7e/wMvnf8jKr7F4olFl48vLmxs+swkpg/+xsUiiXuOaVv3LlQLJGt+ciRff2LQvVUawF6sJltATCzLZIOqrDPBKB8EkdXaOt1l6QS8ADR6Xkrf4yZ9Uj6EDgAeG/gwSXNB+YDHHrooTU+HeecG7Uq5epjqu0zWG6OQ142iWI6SzGVYf/x+/Pq9l0UUhm6UxmmT5nAU+90U8xkKWbaKKYz0b5hfc5JR/KjlS/RncqyM5Xhlu+c8P/FYS6HWutfbAynwKlnMVTvSwLVM/bhHtuLysYYsgCV9ChwSIW7rt7Dv6EKbb3/LZ9vZm9K2o+oAL2QaH7RYI/p32i2FFgK0NnZ6SOkzjm3d/Yk7+5Rbq53Xh63+EZ++JcXKaQyFFJZbrv0hP6FYe+pafWFe9aALxXNHOJv/OO5vv0zR08b6aeQOF7EuZE2ZAFqZqdUu0/SO5I6wuhnB7C1wm5d9J2mh2i+9Zpw7DfD7UeSfkd0Cuie8JhJQJekFPApYNuePCHnnHN7pTfv9poIvFVln6bm5tazvs7ic/qfrs0McXp3OKN3+/oFwF1j7As/l1lPLTU+fgUwN6zPBSpd72AVcKqkcZLGAacCqySlJB0IICkNnAFsrHDcs4HHfP6nc87V1VpgsqTDJWWAOUS5uFwscnO9R+N89M6NhH3h5zLrqdY5oIuAP0iaB7wOnAMgqRO4xMwuNrNtkn5GlNwArgtt7USFaJrokh+P0ne5szuBeyVtJvrvek6NcTrnnBtEmNN5GdGgQSuwzMyelXQdsM7MVuC5ObF8NM7FTU0FqJm9T/9LmvW2rwMuLtteBiwbsM/HRNeSq3TcbkIx65xzrjHM7GHg4QFtPy5b99ycUD4a5+JG+9KZbUnvAq8N4yEHUuGb9QmQ1LghubEnNW5IbuzDjfszZtaUyzW66irl5enTp08vbtiwV8dLT53K+vXr149EbCMsqZ8ziHHsR02delQqnd6rwbKeYrHn6Q0bnh7pmEZIbPt8CCOWl/epAnS4JK0zs85mxzFcSY0bkht7UuOG5Mae1LhdbZL6uic1bkhu7EmNG5Ib+0jGXeuXkJxzzjnnnBsWL0Cdc84551xDjfYCdGmzA9hLSY0bkht7UuOG5Mae1LhdbZL6uic1bkhu7EmNG5Ib+4jFParngDrnnHPOucYb7SOgzjnnnHOuwUZlASpplqRNkjZLWtDseAYjaZKkxyU9L+lZSVeE9vGSHpH0Urgd1+xYK5HUKukpSSvD9uGSngxx/z784krsSBorabmkF0LfH5eEPpd0ZXifbJR0n6S2uPa5pGWStkraWNZWsY8VuSl8ZjdIOrp5kbt68LzcOJ6XG8vzcmWjrgCV1ArcDMwGpgDnSZrS3KgG1QNcZWZHAscCl4Z4FwCrzWwysDpsx9EVwPNl2zcAN4a4PwDmNSWqoS0B/mZmnweOInoOse5zSROAy4FOM/si0a/ZzCG+ff4bYNaAtmp9PBuYHJb5wK0NitE1gOflhvO83CCelwdhZqNqAY4DVpVtLwQWNjuuYcT/EDAT2AR0hLYOYFOzY6sQ68TwZj0JWAmI6AK2qUqvRVwWYH/gFcIc6bL2WPc5MAF4AxhP9CtnK4HT4tznwGHAxqH6GLgdOK/Sfr4kf/G83NBYPS83Nm7Py1WWUTcCSt+boVdXaIs9SYcB04AngYPNbAtAuD2oeZFVtRj4PrArbB8A/MfMesJ2XPv+COBd4K5wmuoOSe3EvM/N7E3gF8DrwBbgQ2A9yejzXtX6OLGfW7dHEvv6el5uGM/LzVOXvDwaC1BVaIv9pQAkjQEeAL5rZv9tdjxDkXQGsNXMyn8yLyl9nwKOBm41s2nAx8TstE4lYV7OmcDhwKeBdqJTJAPFsc+HkpT3jts7iXx9PS83lOfl+KnpvTMaC9AuYFLZ9kTgrSbFskckpYmS3G/N7MHQ/I6kjnB/B7C1WfFVcTzwNUmvAvcTne5ZDIyV1Pu7vnHt+y6gy8yeDNvLiRJf3Pv8FOAVM3vXzIrAg8AMktHnvar1ceI+t25YEvf6el5uOM/LzVOXvDwaC9C1wOTwDbQM0WTgFU2OqSpJAu4EnjezX5bdtQKYG9bnEs1Big0zW2hmE83sMKI+fszMzgceB84Ou8UubgAzext4Q9LnQtPJwHPEvM+JTvEcKykf3je9cce+z8tU6+MVwEXhW5fHAh/2nhJy+wTPyw3gebkpPC9X0+zJrk2aYHs68CLwb+DqZsczRKxfIRrS3gD8KyynE83bWQ28FG7HNzvWQZ7DicDKsH4E8E9gM/BHINvs+KrE/CVgXej3PwPjktDnwLXAC8BG4F4gG9c+B+4jmhNVJPpPel61PiY61XNz+Mw+Q/SN0qY/B19G9P3gebmxz8HzcuPi9rxcYfFfQnLOOeeccw01Gk/BO+ecc865JvIC1DnnnHPONZQXoM4555xzrqG8AHXOOeeccw3lBahzzjnnnGsoL0Cdc84551xDeQHqnHPOOecaygtQ55xzzjnXUP8Dzw3dgfP+o9YAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 792x396 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dm.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tue Jul 23 11:33:23 PDT 2019\r\n"
     ]
    }
   ],
   "source": [
    "!date"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "dismod_mr",
   "language": "python",
   "name": "dismod_mr"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
